Stochastic Vibrations of a System of Plates Immersed in Fluid Using a Coupled Boundary Element, Finite Element, and Finite Difference Methods Approach

The main objective of this work is to investigate the natural vibrations of a system of two thin (Kirchhoff–Love) plates surrounded by liquid in terms of the coupled Stochastic Boundary Element Method (SBEM), Stochastic Finite Element Method (SFEM), and Stochastic Finite Difference Method (SFDM) implemented using three different probabilistic approaches. The BEM, FEM, and FDM were used equally to describe plate deformation, and the BEM was applied to describe the dynamic interaction of water on a plate surface. The plate’s inertial forces were expressed using a diagonal or consistent mass matrix. The inertial forces of water were described using the mass matrix, which was fully populated and derived using the theory of double-layer potential. The main aspect of this work is the simultaneous application of the BEM, FEM, and FDM to describe and model the problem of natural vibrations in a coupled problem in solid–liquid mechanics. The second very important novelty of this work is the application of the Bhattacharyya relative entropy apparatus to test the safety of such a system in terms of potential resonance. The presented concept is part of a solution to engineering problems in the field of structure and fluid dynamics and can also be successfully applied to the analysis of the dynamics of the control surfaces of ships or aircraft.


Introduction
One of the computational numerical tools is the Boundary Element Method (BEM), which may be applied in the numerical analysis of a variety of engineering and scientific problems, e.g., thin or thick plate bending. The mathematical foundation of this method began with the theory of potential developed by Kellog [1]. Following, Kupradze [2] and Michlin [3] proposed multidimensional and singular boundary integral equations for selected parts of the theory of elasticity. Further works were carried out by Jaswon [4], Symm [5], Jaswon, and Ponter [6], who proposed a direct application of the integral equations to static potential problems. The above investigation can be treated as the theoretical foundations of the boundary integral equation method. The BEM is an independent tool that may be optimal for selected engineering problems. The main advantage of the BEM approach is its simplicity in formulating a large number of structural mechanics problems based on the theories of potential or elasticity. For some problems, this method is irreplaceable. A significant degree of discretization is not required using the BEM approach. The problem of plate bending was solved using the BEM by Stern [7], Altiero and Sikarskie [8], and Bèzine and Gamby [9]. The effectiveness of the BEM and FEM was compared by Debbich [10,11] for the stress analysis of plates of medium thickness. A modification of the Kirchhoff-Love plate theory in terms of the BEM was performed by El-Zafrany, Debbih,

Fluid-Structure Inertial Interaction for a Single Plate
It was assumed that the considered fluid is incompressible and inviscid. There was no separate flow between the plate and fluid, and no fluid moved along the plane direction of the plate. The surrounding fluid is the source of the additional inertial forces.
It was assumed that the velocity potential of a fluid for small disturbances is expressed as where t represents time, x = (x, y, z) is the position vector of a point P arbitrarily assumed in space, and ω expresses the circular frequency of a fluid. It was assumed that the velocity potential of liquid (1) must satisfy the Laplace equation.
The solution of Equation (2) can be expressed as the double-layer potential in terms of the boundary integral equation. wherein is the fundamental solution of the Laplace Equation (2), i = √ −1, ∼ V 1 (Q) and ∼ V 2 (Q) are the amplitudes of the velocity potential above and below the surface, and ∼ V(P) is the amplitude of the velocity potential in any point of space. Then, it was assumed that point P approaches surface S ( Figure 1).
Subsequently, the velocity potential first derivative (3) specified in the direction normal to surface S at any point P is calculated Subsequently, the velocity potential first derivative (3) specified in the direction normal to surface S at any point P is calculated The considered Neumann's boundary conditions allows connection between the fluid velocity potential and the displacement of surface S by a differential relation wherein ( , ) = ( ) is a displacement that is perpendicular to the plate domain. Then, the hydrodynamic pressure acting on plate surface S can be determined and given by the fluid's density coupled with the first derivative of the fluid velocity potential with respect to time Now, Equation (5) can be differentiated again with respect to time, and the above expressions (6) and (7) are taking place for it. The above operations lead to the relation between the circular frequencies of a fluid and amplitude of displacement with the amplitude of the resultant hydrodynamic pressure at point Q placed on surface S − (P) = Δ (Q) * (P,Q) , wherein (Q) = (Q) − (Q) is the difference in hydrodynamic pressure between the upper and lower surfaces of the plate. The boundary integral in Equation (8) can be solved using the Boundary Element Method (BEM). The fluid surface adjacent to the plate is divided into a finite number of subsurfaces that act as boundary elements for the fluid. Using the simplest approximation for a rectangular plate, each subsurface has only one collocation point, which is placed centrally (Figure 2). The plate domain is divided into a finite number of subdomain surfaces with areas Sn, and the amplitude of displacement in arbitrary points (xm, ym) are coupled with the hydrodynamic pressure amplitude The considered Neumann's boundary conditions allows connection between the fluid velocity potential and the displacement of surface S by a differential relation wherein w(x, t) = ∼ w(x)e iωt is a displacement that is perpendicular to the plate domain. Then, the hydrodynamic pressure acting on plate surface S can be determined and given by the fluid's density coupled with the first derivative of the fluid velocity potential with respect to time Now, Equation (5) can be differentiated again with respect to time, and the above expressions (6) and (7) are taking place for it. The above operations lead to the relation between the circular frequencies of a fluid and amplitude of displacement with the amplitude of the resultant hydrodynamic pressure at point Q placed on surface S −ρ f ω 2 ∼ w(P) = S ∆ ∼ p(Q) wherein ∆ ∼ p(Q) = ∼ p 2 (Q) − ∼ p 1 (Q) is the difference in hydrodynamic pressure between the upper and lower surfaces of the plate.
The boundary integral in Equation (8) can be solved using the Boundary Element Method (BEM). The fluid surface adjacent to the plate is divided into a finite number of subsurfaces that act as boundary elements for the fluid. Using the simplest approximation for a rectangular plate, each subsurface has only one collocation point, which is placed centrally ( Figure 2). Subsequently, the velocity potential first derivative (3) specified in the direction normal to surface S at any point P is calculated The considered Neumann's boundary conditions allows connection between the fluid velocity potential and the displacement of surface S by a differential relation wherein ( , ) = ( ) is a displacement that is perpendicular to the plate domain. Then, the hydrodynamic pressure acting on plate surface S can be determined and given by the fluid's density coupled with the first derivative of the fluid velocity potential with respect to time Now, Equation (5) can be differentiated again with respect to time, and the above expressions (6) and (7) are taking place for it. The above operations lead to the relation between the circular frequencies of a fluid and amplitude of displacement with the amplitude of the resultant hydrodynamic pressure at point Q placed on surface S − (P) = Δ (Q) * (P,Q) , wherein (Q) = (Q) − (Q) is the difference in hydrodynamic pressure between the upper and lower surfaces of the plate. The boundary integral in Equation (8) can be solved using the Boundary Element Method (BEM). The fluid surface adjacent to the plate is divided into a finite number of subsurfaces that act as boundary elements for the fluid. Using the simplest approximation for a rectangular plate, each subsurface has only one collocation point, which is placed centrally ( Figure 2). The plate domain is divided into a finite number of subdomain surfaces with areas Sn, and the amplitude of displacement in arbitrary points (xm, ym) are coupled with the hydrodynamic pressure amplitude The plate domain is divided into a finite number of subdomain surfaces with areas S n , and the amplitude of displacement in arbitrary points (x m , y m ) are coupled with the hydrodynamic pressure amplitude The boundary integral Equation (9) can be rewritten in the matrix notation to the form −4πρ f ω 2 ∼ w = H ∼ p (10) and H is the (N × N) matrix with elements derived as the integral over the nth subsurface, wherein N expresses the total number of internal subsurfaces coupled with single collocation points (subsurfaces of the constant type) [22,23].
Now, the hydrodynamic forces acting on the plate surface can be represented as: where P is the vector of the resultant hydrodynamic forces, where each element is defined as P n = ∆p n S n , and M f is the fluid mass matrix. This matrix is expressed by the relation with the diagonal matrix S which groups areas of internal subdomains All of the elements of the matrix H can be calculated analytically, e.g., [22,46], according to designations presented in Figure 2 H mn = − 1 Next, a description of the problem of two plates surrounded by liquid is presented. The problem of the dynamics of a single plate interacting with a liquid can be extended to the problem of the dynamics of a finite number of plates immersed in a liquid medium.
A medium constraint in the form of a rigid partition, tank bottom, or free surface may also be assigned to such a system. This issue is illustrated using an example of natural and forced vibrations of a system of two plates immersed in a liquid medium. A system of two plates immersed in a liquid medium is presented. The plates were situated in the same plane ( Figure 3) or one above another ( Figure 4).
The boundary integral Equation (9) can be rewritten in the matrix notation to the form and H is the (N × N) matrix with elements derived as the integral over the nth subsurface, wherein N expresses the total number of internal subsurfaces coupled with single collocation points (subsurfaces of the constant type) [22,23].
Now, the hydrodynamic forces acting on the plate surface can be represented as: where P is the vector of the resultant hydrodynamic forces, where each element is defined as = , and Mf is the fluid mass matrix. This matrix is expressed by the relation with the diagonal matrix S which groups areas of internal subdomains All of the elements of the matrix H can be calculated analytically, e.g., [22,46], according to designations presented in Figure Next, a description of the problem of two plates surrounded by liquid is presented. The problem of the dynamics of a single plate interacting with a liquid can be extended to the problem of the dynamics of a finite number of plates immersed in a liquid medium.
A medium constraint in the form of a rigid partition, tank bottom, or free surface may also be assigned to such a system. This issue is illustrated using an example of natural and forced vibrations of a system of two plates immersed in a liquid medium. A system of two plates immersed in a liquid medium is presented. The plates were situated in the same plane ( Figure 3) or one above another ( Figure 4).   Therefore, it is necessary to construct the matrix H in the modified form: where matrices (P,Q) , (P,Q ') , (P ',Q) , and (P ',Q ′) consist of surface integrals from the fundamental solution (12) over the subsurfaces belonging to the first or the second plate.
In view of the above, the fluid mass matrix for the system of two plates interacting with a liquid is defined: where ( , ) groups the areas of the interior subsurfaces of the first and second plates, respectively.

The Boundary Element Method Solution
A thin (Kirchhoff-Love) plate is surrounded by a fluid, and its vibrations have small amplitudes. Inside the plate domain, additional collocation points, coupled with lumped masses, were introduced according to Bèzine's approach [9]. In each collocation occurring inside a plate domain, vectors of displacements ( ) , accelerations ( ) , and inertial forces ( ) dependent on time t were introduced, e.g., [22,23] ( ) = , where the amplitude of a plate's inertial force is given by = , is a lumped mass, is the amplitude of displacement, and is the natural frequency of the plate vibrations.
The boundary and domain integral equations describing the plate deformation have the characteristics of amplitude equations in terms of the modified approach of a thin plate bending description, e.g., [45]. These equations may be derived using, e.g., Betti's theorem, or otherwise they are integral representations of the thin plate differential equation solution Therefore, it is necessary to construct the matrix H in the modified form: where matrices H(P, Q), H P, Q , H P , Q , and H P , Q consist of surface integrals from the fundamental solution (12) over the subsurfaces belonging to the first or the second plate.
In view of the above, the fluid mass matrix for the system of two plates interacting with a liquid is defined: where S (1,2) groups the areas of the interior subsurfaces of the first and second plates, respectively.

The Boundary Element Method Solution
A thin (Kirchhoff-Love) plate is surrounded by a fluid, and its vibrations have small amplitudes. Inside the plate domain, additional collocation points, coupled with lumped masses, were introduced according to Bèzine's approach [9]. In each collocation occurring inside a plate domain, vectors of displacements w i (t), accelerations .. w i (t), and inertial forces B i (t) dependent on time t were introduced, e.g., [22,23] where the amplitude of a plate's inertial force is given by ∼ B i = ω 2 m i W i , m i is a lumped mass, W i is the amplitude of displacement, and ω is the natural frequency of the plate vibrations.
The boundary and domain integral equations describing the plate deformation have the characteristics of amplitude equations in terms of the modified approach of a thin plate bending description, e.g., [45]. These equations may be derived using, e.g., Betti's theorem, where w * (y, x) = (1/(8πD))r 2 ln(r) is the Green function, which is the fundamental solution of the biharmonic equation ∇ 4 w * (y, x) = (1/D)δ(y, x) for a thin isotropic plate, δ is the Dirac delta, r = |y − x|, x is the source point, y is a field point, D = Eh 3 / 12 1 − v 2 is the plate stiffness, h is the plate thickness, and E and v are Young's modulus and Poisson's ratio, respectively. The relation between the angle of rotation φ s (y) and deflection w(y) is specified by the differential relationship φ s (y) = dw(y)/ds, e.g., [22]. The Kirchhoff corner forces were expressed (replaced) by the forces distributed along the boundary elements  (19) and (20) is responsible for the coupled inertial forces of a plate-fluid system, e.g., [22,23,45]. The value of the coefficient c(x) occurring in Equations (20) and (21) is 1 for x located inside the plate domain, 0.5 for x located on the smooth boundary, and 0 for x located outside the plate domain.
The BEM methodology was applied with the discretization of a plate edge using boundary elements, in this case, of the constant type, and a plate domain was divided into rectangular subsurfaces coupled with one collocation point as shown in Figure 5. Each of these plate subsurfaces was coupled with the corresponding ones describing the acting liquid.
where * ( , ) = (1 (8 ) ⁄ ) ( ) is the Green function, which is the fundamental solution of the biharmonic equation ∇ * ( , ) = (1⁄ ) ( , ) for a thin isotropic plate, is the Dirac delta, = | − | , x is the source point, y is a field point, = ℎ 12(1 − ) ⁄ is the plate stiffness, h is the plate thickness, and E and v are Young's modulus and Poisson's ratio, respectively. The relation between the angle of rotation ( ) and deflection ( ) is specified by the differential relationship ( ) = ( )⁄ , e.g., [22]. The Kirchhoff corner forces were expressed (replaced) by the forces distributed along the boundary elements placed in the vicinity of the plate corners. Variables * ( , ), * ( , ) , * ( , ) , * ( , ) , * ( , ) , and * ( , ) denote the fundamental functions were derived directly from the fundamental solution * ( , ) as its derivatives, according to the thin plate (Kirchhoff-Love) theory. The jth amplitude of the resultant inertial forces = + occurring in Equations (19) and (20) is responsible for the coupled inertial forces of a plate-fluid system, e.g., [22,23,45]. The value of the coefficient ( ) occurring in Equations (20) and (21) is 1 for x located inside the plate domain, 0.5 for x located on the smooth boundary, and 0 for x located outside the plate domain.
The BEM methodology was applied with the discretization of a plate edge using boundary elements, in this case, of the constant type, and a plate domain was divided into rectangular subsurfaces coupled with one collocation point as shown in Figure 5. Each of these plate subsurfaces was coupled with the corresponding ones describing the acting liquid.
groups the amplitudes of the boundary quantities: The vector is formed by the amplitudes of the internal displacements of the collocation points and where matrices wB contain appropriate boundary integrals for a single plate (j = 1, 2) depending on the type of boundary. In the matrix G and M (j) N , λ = ω 2 , I is the unit matrix, and N is the number of lumped masses for a single plate. Elimination of the boundary variables ∼ X (1,2) from the matrix Equation (22) leads to the standard eigenvalue problem where ∼ λ = 1/ω 2 and Note the fundamental function used to compute the matrix H (1,2) elements. This function is obtained by performing a double differentiation of the function 1/r with respect to the z coordinate, where r = x 2 + y 2 + z 2 [22,45].

The Finite Element Method Approach to the Fluid-Structure Interaction
The free vibration problem of the structure can be described by a generalized eigenvalue problem, which can be written in matrix notation as follows where K (1,2) and M (1,2) are the stiffness and mass matrices of the whole system, respectively, ω is the eigenvalue, and ∼ w (1,2) is the nonzero eigenvector. The matrices K (1,2) and M (1,2) have the form where K (1) and K (2) are the stiffness matrices of the first and second plates, respectively, and the mass matrix of the plate-fluid system is coupled with the fluid according to relation (17).
Vector ∼ w (1,2) has already been defined in Equation (24). The bending of a single plate with constant thickness is described by a rectangular four-node finite element with three degrees of freedom at each node. At each element i th node, in the Cartesian coordinate system, there are introduced: deflection w i and two angles of rotation in mutually perpendicular directions ϕ ix and ϕ iy , respectively, wherein the function of deflection is expressed as the polynomial of the fourth order [51] w(x, y) = α 1 where α 11 = 0 or α 12 = 0. The displacement field inside the finite element is presented in the FEM methodology i expresses the number of the actual node, and j expresses the current degree of freedom at the ith node. The set of shape functions N (i) j for each node can be given by formulas derived directly from relation (32). A detailed description of the considered finite element was presented by Kuczma [51] and quoted by Lenartowicz and Guminiak [48].
The stiffness matrix of the finite element is defined according to FEM methodology where h is the plate thickness, B is the shape functions derivatives matrix, and D is the plate physical quantities matrix. The mass matrix of the single element M e has the characteristics of the consistent matrix. The boundary conditions for the plates were introduced according to the well-known FEM methodology.

The Finite Difference Method Solution to the Fluid-Structure Interaction
The free vibration in terms of the Finite Difference Method is described by the difference procedures leading to the following matrix equation which has the characteristics of the generalized eigenvalue problem The matrix of suitable difference operators for the system of two plates K and group matrices of the difference operators of the first and second plates K (1) ∆ and K (2) ∆ ( Figure 6) The mass matrix of the plate-fluid system was coupled with the fluid according to relation (17); the mass matrix of a single plate M is a diagonal matrix containing plate masses m i , which are concentrated in subsequent FDM mesh nodes, divided by their surface area. Finally, the vector ∼ w (1,2) is expressed by the relation (24).

= ( )
and group matrices of the difference operators of the first and second plates ( ) and ( ) ( Figure 6) The mass matrix of the plate-fluid system was coupled with the fluid according to relation (17); the mass matrix of a single plate M is a diagonal matrix containing plate masses mi, which are concentrated in subsequent FDM mesh nodes, divided by their surface area. Finally, the vector ( , ) is expressed by the relation (24).  The boundary conditions for the plates are introduced according to the FDM methodology, described comprehensively [48].

Probabilistic Responses and Relative Entropy Application
The Gaussian uncertainty in some of the materials or geometrical design parameters of the aforementioned plate was assumed. Random quantities include the thickness of the plate, Young's modulus, Poisson's ratio, distance of the collocation point from the plate edge, boundary conditions, and others. The structural response functions were realized using polynomial curves of the third order. Structural uncertainty numerical analysis of the considered plates was performed using the semi-analytical method (SAM), iterative generalized 10th order stochastic perturbation technique (SPT), and Monte Carlo simulation (MCS).
First, for the SAM approach, the probability distribution of the observed quantity must be assumed in advance. This can be, for example, a Gaussian distribution. Next, the polynomial expressions should be applied directly to the probabilistic integral formulas, which allows the analytical calculation of all of the formulas for the expected values, coefficients of variation, etc., to have been implemented due to the derivations presented in [49]. This approach is relatively simple and easy to implement.
A more complex approach is the Stochastic Perturbation Technique (SPT). In this study, the Taylor series expansion of all the input variables and state functions was used, assuming that the perturbation parameter is ε. This expansion is then used to determine random moments. Random moments can be calculated analytically or determined in symbolic form if and only if the function of the random parameter is known.
The Monte Carlo simulation approach appears to be the simplest because it is based on the so-called almost large numbers; unfortunately, it entails a large number of operations (trials) on the order of fifty or one hundred thousand. However, in the case of a classical numerical simulation, a certain probability distribution of the occurrence of a random quantity should be assumed. Therefore, this method is characterized by a high time cost.
R denotes the admissible limit of the given structure, and E denotes its extreme effort. The previous engineering design codes allowed us to make the following interpretations in the case of the eigenfrequency and extreme excitation. The distance in between them cannot be smaller than a quarter of this eigenfrequency. This was carried out to avoid structural resonance. Therefore, the satisfactory safety of the given system may be measured using the following FORM reliability index: where ω i stands for each next eigenfrequency. Quite a similar calculation can be provided with the use of a relative entropy H. This entropy can be calculated due to the Bhattacharyya theory as This can be rewritten as The final safety measure compensable to β form is obtained as

Numerical Examples
The natural vibrations of the system of two thin elastic plates surrounded by liquid on all sides were considered. The system of plates was analyzed numerically using three independent Stochastic Boundary Element Methods (SBEM). Numerical comparison of the deterministic BEM, FEM, and FDM solutions for the first four fundamental natural frequencies for a plate immersed in fluid were performed as the foundations for further SBEM, SFEM, and SFDM procedures with Young's modulus equal to E = 205 GPa, Poisson's ratio v = 0.3, plate dimensions l x = l y = 2.0 m, plate thickness h = 0.05 m, and fluid density ρ f = 1000 kg/m 3 . The numerical results were compared with the analytical solutions obtained by C.C. Liang et al. [44] and Lenartowicz and Guminiak [48] for a single plate fully submerged in fluid.
The BEM-BEM numerical computations were carried out using the original computational program developed by the first author of the Fortran 90 programming language. The direct version of the BEM and the static fundamental solution of the Kirchhoff plate problem were used, where the quasi-diagonal boundary terms in the characteristic matrix were calculated analytically, and the nondiagonal terms were calculated numerically by applying a 12-point Gauss quadrature. The angle of rotation in the tangent direction occurring on the free edge was calculated by formulating the difference relations between displacements for three boundary nodes (collocation points) taking place next to each other.
The FEM-BEM and FDM-BEM deterministic computations were processed using the numerical procedures developed by the last two authors of the OCTAVE numerical package. Finally, probabilistic analyses were performed using the original procedures developed by the second author of the MAPLE v.19 scientific package, where LSM fittings and three probabilistic approaches (SAM, SPT, and MCS) were implemented. The following basic discretization of a single plate was taken: the number of boundary elements was 100, the number of internal subdomains was 100, the number of finite elements was 256, and the number of finite difference grids was 900 for the BEM-BEM, FEM-BEM, and FDM-BEM approaches, respectively.
For probabilistic investigation, the following intervals of variables and their representations for the needs of the LSM analysis are presented:  Table 1. It is clear from these results that the BEM-BEM approach generally returns all the basic natural frequencies larger than the analytical method, yet still smaller than the combined FEM-BEM and FDM-BEM approaches. Moreover, an influence of the parameter ∼ δ is rather marginal because even dramatic changes of its mean value result in negligible modifications of the aforementioned frequencies. This is an important result, which shows no sensitivity from the BEM approach towards the numerical parameter of the plate(s) discretization.

Two Cantilever Plates Immersed in a Fluid and Placed in One Plane along One Line
First, the results of two modes of natural vibrations of the system of cantilever plates fully submerged in water from the BEM-BEM solution and a = 1.0 m are shown in Figure 7.
The influence of parameter ∼ δ on the structural response for the first four natural frequencies of the two cantilever plates fully submerged in fluid for the BEM-BEM solution is presented in Table 2. A comparison of the results from the FEM-BEM and FDM-BEM solutions is presented in Table 3 for basic discretization using a 16

Two Cantilever Plates Immersed in a Fluid and Placed in One Plane along One Line
First, the results of two modes of natural vibrations of the system of cantilever plates fully submerged in water from the BEM-BEM solution and a = 1.0 m are shown in Figure  7. The influence of parameter on the structural response for the first four natural frequencies of the two cantilever plates fully submerged in fluid for the BEM-BEM solution is presented in Table 2. A comparison of the results from the FEM-BEM and FDM-BEM solutions is presented in Table 3 Table 3. Natural frequencies for the system of two cantilever plates fully immersed in water and a = 1.0 m obtained for the FEM-BEM and FDM-BEM approaches.  The results presented in these tables lead to almost the same conclusion as the data included in Table 1; not a single plate, nor the system of two plates in fluid show any sensitivity to the parameter

Stochastic BEM-BEM Analysis
The main aim of the first numerical experiment was to compare the proposed triple probabilistic methods: the semi-analytical method (SAM), perturbative technique (SPT), and Monte Carlo simulations (MCS). The set of deterministic problems related to the system of two plates immersed in fluid was solved using the BEM to describe the plate bending and to determine the inertial forces coming from the fluid (BEM-BEM analysis). A comparison was performed for the random placement of the distance between these plates expressed by parameter a for a plate exhibiting a Gaussian uncertainty natural frequency defined by its mean value and the specific range of its coefficient of variation, i.e., α(b) ∈ [0.00, 0.025]. The response function for the first natural frequency is expressed as the third-order polynomial: The response function for the second natural frequency is expressed as the third-order polynomial: The results of the calculation of the expected values E(ω), coefficients of variation α(ω), skewness β(ω), and kurtosis κ(ω) of the first and the second natural frequencies are presented in Figure 8. First of all, it was observed that the three probabilistic methods return similar results each time; therefore, they can be used alternatively as a solution to this problem. Further, the influence of the uncertainty within this parameter has a negligibly small impact on the statistical scattering of the natural frequencies of the given system; the resulting coefficients of variation equal all almost 0. Finally, the resulting eigenfrequency probability distribution is rather distant from the Gaussian one because higher-order statistics clearly diverge from 0 while increasing the coefficient CoV(a).
Further, the polynomial response functions for the first two natural frequencies corresponding to varying fluid density are expressed by the third-order polynomials: The results of the basic probabilistic characteristics calculation, that is, the expected values ( ) , coefficients of variation ( ) , skewness ( ) , and kurtosis ( ) of the first and the second natural frequencies, are presented in Figure 9 as functions of the input coefficient of variation in fluid density. This situation is remarkably different than in the First of all, it was observed that the three probabilistic methods return similar results each time; therefore, they can be used alternatively as a solution to this problem. Further, the influence of the uncertainty within this parameter has a negligibly small impact on the statistical scattering of the natural frequencies of the given system; the resulting coefficients of variation equal all almost 0. Finally, the resulting eigenfrequency probability distribution is rather distant from the Gaussian one because higher-order statistics clearly diverge from 0 while increasing the coefficient CoV(a).
Further, the polynomial response functions for the first two natural frequencies corresponding to varying fluid density are expressed by the third-order polynomials: The results of the basic probabilistic characteristics calculation, that is, the expected values E(ω), coefficients of variation α(ω), skewness β(ω), and kurtosis κ(ω) of the first and the second natural frequencies, are presented in Figure 9 as functions of the input coefficient of variation in fluid density. This situation is remarkably different than in the case of plate thickness randomness; although the first two moments are equal when computed using different probabilistic strategies, a higher order statistics agreement is obtained for the first eigenfrequency only. Fluid density has a tremendous impact on natural frequencies, and their coefficients of variation are remarkably larger than the coefficient of variation of the fluid density.
The results of the calculation with the expected values ( ), coefficients of variation ( ), skewness ( ), and kurtosis ( ) of the first and second natural frequencies are presented in Figure 10. This parameter of uncertainty had almost no statistical influence on the natural frequencies, and the PDFs of these frequencies were remarkably different from the Gaussian one.  The results of the calculation with the expected values E(ω), coefficients of variation α(ω), skewness β(ω), and kurtosis κ(ω) of the first and second natural frequencies are presented in Figure 10. This parameter of uncertainty had almost no statistical influence on the natural frequencies, and the PDFs of these frequencies were remarkably different from the Gaussian one. The next calculations were performed using a random approach on the restrained support of both plates. The propagation of the unsupported segment is described by the parameter "s", as shown in Figure 11. The next calculations were performed using a random approach on the restrained support of both plates. The propagation of the unsupported segment is described by the parameter "s", as shown in Figure 11. The next calculations were performed using a random approach on the restrained support of both plates. The propagation of the unsupported segment is described by the parameter "s", as shown in Figure 11.  Figure 11. The propagation of the unsupported segment is described by the parameter "s".
The response functions of the first two natural frequencies for the parameter "s" describing the propagation of the unsupported section of the plates are expressed by the third-order polynomials: ω 2 (s) = 38.5197594349883 − 4.76747491350958 · s− −18.2226428613072 · s 2 − 10.5629477770474·s 3 .
The results of the calculation of the expected values E(ω), coefficients of variation α(ω), skewness β(ω), and kurtosis κ(ω) of the first and second natural frequencies are presented in Figure 12. The results of the calculation of the expected values ( ), coefficients of variation ( ), skewness ( ), and kurtosis ( ) of the first and second natural frequencies are presented in Figure 12. Now, once more SAM, SPT, and MCS probabilistic methods coincide when all characteristics are considered, whereas the uncertainty in the plate separation parameter s has a rather marginal influence on the random output characteristics. This impact is larger than in the case of the parameter , but definitely smaller than in the case of the fluid density.
Finally, the plate thickness uncertainty was studied, where the response functions of the first two natural frequencies were approximated via the Least Squares Method by the third-order polynomials: The results of the calculation with the expected values ( ), coefficients of variation ( ), skewness ( ), and kurtosis ( ) of the first and second natural frequencies are presented in Figure 13. There is no doubt that this parameter should play, and really does play, a very important role in the uncertainty analysis because the output CoV is somewhat larger than that defined for this thickness. A coincidence of the MCS statistical estimation is discussive in case of the kurtosis only; nevertheless, the resulting eigenfrequency PDFs are close to the Gaussian density parameters. Now, once more SAM, SPT, and MCS probabilistic methods coincide when all characteristics are considered, whereas the uncertainty in the plate separation parameter s has a rather marginal influence on the random output characteristics. This impact is larger than in the case of the parameter ∼ δ, but definitely smaller than in the case of the fluid density. Finally, the plate thickness uncertainty was studied, where the response functions of the first two natural frequencies were approximated via the Least Squares Method by the third-order polynomials: The results of the calculation with the expected values E(ω), coefficients of variation α(ω), skewness β(ω), and kurtosis κ(ω) of the first and second natural frequencies are presented in Figure 13. There is no doubt that this parameter should play, and really does play, a very important role in the uncertainty analysis because the output CoV is somewhat larger than that defined for this thickness. A coincidence of the MCS statistical estimation is discussive in case of the kurtosis only; nevertheless, the resulting eigenfrequency PDFs are close to the Gaussian density parameters. The results of the calculation with the expected values ( ), coefficients of variation ( ), skewness ( ), and kurtosis ( ) of the first and second natural frequencies are presented in Figure 14.

Stochastic FEM-BEM Analysis
In this study, the main aim of the second numerical experiment was to check probabilistic approaches by comparing the SAM, SPT, and MCS results. The deterministic problem of the system of two plates immersed in fluid was solved by using the FEM to describe plate bending and the BEM to describe the occurrence of fluid (FEM-BEM analysis). The response functions of the first two natural frequencies are expressed by the third-order polynomials: The results of the calculation with the expected values E(ω), coefficients of variation α(ω), skewness β(ω), and kurtosis κ(ω) of the first and second natural frequencies are presented in Figure 14.
The three probabilistic methods coincided with each other for all probabilistic characteristics, and the induced uncertainty increased, almost linearly, while the input coefficient variation increased. The remaining characteristic of the structural response demonstrated a very similar, yet nonlinear, monotonous increase. It should be mentioned that the impact of the fluid density's uncertainty is smaller than the one verified above, implying that the similarity of the mean values in BEM-BEM and FEM-BEM approaches does not imply any similarity in higher-order statistics. Finally, it should be pointed out that these results are extremely similar to the ones presented in Figure 14, which were obtained while replacing the FEM with the FDM. The three probabilistic methods coincided with each other for all probabilistic characteristics, and the induced uncertainty increased, almost linearly, while the input coefficient variation increased. The remaining characteristic of the structural response demonstrated a very similar, yet nonlinear, monotonous increase. It should be mentioned that the impact of the fluid density's uncertainty is smaller than the one verified above, implying that the similarity of the mean values in BEM-BEM and FEM-BEM approaches does not imply any similarity in higher-order statistics. Finally, it should be pointed out that these results are extremely similar to the ones presented in Figure 14, which were obtained while replacing the FEM with the FDM.

Stochastic FDM-BEM Analysis
Similarly, the main aim of the third numerical experiment was to validate the proposed triple probabilistic methodology by comparing the SAM, SPT, and MC approaches. The deterministic problem of the system of two plates immersed in fluid was solved using the FDM description of the plate bending and the BEM description of the occurrence fluid (FDM-BEM analysis). The response functions of the first two natural frequencies are expressed by the third-order polynomials: The results of the calculation with the expected values ( ), coefficients of variation ( ), skewness ( ), and kurtosis ( ) of the first and the second natural frequencies are presented in Figure 15.

Stochastic FDM-BEM Analysis
Similarly, the main aim of the third numerical experiment was to validate the proposed triple probabilistic methodology by comparing the SAM, SPT, and MC approaches. The deterministic problem of the system of two plates immersed in fluid was solved using the FDM description of the plate bending and the BEM description of the occurrence fluid (FDM-BEM analysis). The response functions of the first two natural frequencies are expressed by the third-order polynomials: The results of the calculation with the expected values E(ω), coefficients of variation α(ω), skewness β(ω), and kurtosis κ(ω) of the first and the second natural frequencies are presented in Figure 15.

Two Cantilever Plates Immersed in Fluid and Placed One Directly above the Other
Two modes of natural vibrations of a system of cantilever plates fully submerged in fluid (water) for the BEM-BEM solution and c = 1.0 m are shown in Figure 16.
The influence of the parameter ∼ δ on the structural response for the first four natural frequencies of the two cantilever plates surrounded by water for the BEM-BEM solution is presented in Table 6. Further, a comparison of the results for the FEM-BEM and FDM-BEM solutions and the basic discretization using a 16 × 16 plate FEM mesh with a 15 × 15 fluid BEM mesh and a 30 × 30 plate FDM grid with a 15 × 15 fluid BEM mesh are presented in Table 7. The results showed almost the same properties as in the previous section, an insensitivity to the parameter ∼ δ, whereas an increase in the value of the natural frequencies were obtained in turn for the BEM-BEM, FEM-BEM as well as for the FDM-BEM coupled numerical strategies.

Two Cantilever Plates Immersed in Fluid and Placed One Directly above the Other
Two modes of natural vibrations of a system of cantilever plates fully submerged in fluid (water) for the BEM-BEM solution and c = 1.0 m are shown in Figure 16. The influence of the parameter on the structural response for the first four natural frequencies of the two cantilever plates surrounded by water for the BEM-BEM solution is presented in Table 6.  Table 7. The results showed almost the same properties as in the previous section, an insensitivity to the parameter , whereas an increase in the value of the natural frequencies were obtained in turn for the BEM-BEM, FEM-BEM as well as for the FDM-BEM coupled numerical strategies.

Two Cantilever Plates Immersed in Fluid and Placed One Directly above the Other
Two modes of natural vibrations of a system of cantilever plates fully submerged in fluid (water) for the BEM-BEM solution and c = 1.0 m are shown in Figure 16. The influence of the parameter on the structural response for the first four natural frequencies of the two cantilever plates surrounded by water for the BEM-BEM solution is presented in Table 6 Table 7. The results showed almost the same properties as in the previous section, an insensitivity to the parameter , whereas an increase in the value of the natural frequencies were obtained in turn for the BEM-BEM, FEM-BEM as well as for the FDM-BEM coupled numerical strategies. Figure 16. Two modes of two cantilever plates immersed in fluid and placed one above another for BEM-BEM solution. Table 6. An influence of the parameter ∼ δ on the natural frequency ω for the system of two cantilever plates completely immersed in water and a = 1.0 m.  Table 7. Natural frequencies ω for the system of two cantilever plates fully immersed in water and c = 1.0 m obtained for FEM-BEM and FDM-BEM approaches.  Tables 8 and 9 show comparisons of the results obtained from three different plate FEM meshes, plate FDM grids, and fluid BEM meshes. First, a comparison has been performed for the random distance placement c between two plates having Gaussian distribution uniquely defined by its mean value, and a specific range of the coefficient of variation, i.e., [0.00, 0.025]. The response functions of the first two natural frequencies were obtained as the third-order polynomials:

FEM-BEM-Present FDM-BEM-Present
The expected values E(ω), coefficients of variation α(ω), skewness β(ω), and kurtosis κ(ω) of the first and second natural frequencies are presented in Figure 17. A coincidence of the given three probabilistic numerical techniques is noticeable. Both obtained frequencies have apparent non-Gaussian distributions, and their statistical dispersion is rather marginal in this particular case. coincidence of the given three probabilistic numerical techniques is noticeable. Both obtained frequencies have apparent non-Gaussian distributions, and their statistical dispersion is rather marginal in this particular case. The results of the calculation with the expected values ( ), coefficients of variation ( ), skewness ( ), and kurtosis ( ) of the first and second natural frequencies are presented in Figure 18. The resulting uncertainty is a little bit more influential, and these probabilistic results agree very well with those presented in Section 5.1. Probabilistic damping represents the ratio of output to input uncertainty and it equals about 3.0 here. The next calculations were conducted using a random fluid density, where the response functions of the first two natural frequencies were obtained in the form of the third-order polynomials: The results of the calculation with the expected values E(ω), coefficients of variation α(ω), skewness β(ω), and kurtosis κ(ω) of the first and second natural frequencies are presented in Figure 18. The resulting uncertainty is a little bit more influential, and these probabilistic results agree very well with those presented in Section 5.1. Probabilistic damping represents the ratio of output to input uncertainty and it equals about 3.0 here.

Stochastic FEM-BEM Analysis
Similarly, a probabilistic analysis was performed using the random placement of the distance between the plates expressed by parameter c, the plate exhibiting a Gaussian uncertainty natural frequency by its mean value, and a specific range of its coefficient of variation, i.e., α(b) ∈ [0.00, 0.025]. The response functions of the first two natural frequencies were obtained as the third-order polynomials:

Stochastic FEM-BEM Analysis
Similarly, a probabilistic analysis was performed using the random placement of the distance between the plates expressed by parameter c, the plate exhibiting a Gaussian uncertainty natural frequency by its mean value, and a specific range of its coefficient of variation, i.e., ( ) ∈ 0.00, 0.025 . The response functions of the first two natural frequencies were obtained as the third-order polynomials: The results of the calculation with the expected values ( ), coefficients of variation ( ), skewness ( ), and kurtosis ( ) of the first and the second natural frequencies are presented in Figure 19. Upon further reliability assessment, this parameter appeared to be negligible, whereas the probabilistic results coincided almost perfectly with those obtained in Section 5.2.1. The results of the calculation with the expected values E(ω), coefficients of variation α(ω), skewness β(ω), and kurtosis κ(ω) of the first and the second natural frequencies are presented in Figure 19. Upon further reliability assessment, this parameter appeared to be negligible, whereas the probabilistic results coincided almost perfectly with those obtained in Section 5.2.1.

Stochastic FDM-BEM Analysis
Similarly, the basic purpose of the third numerical experiment was to validate the proposed triple probabilistic methodology: the SAM, SPT, and MC approaches. The deterministic problem of the system of two plates immersed in fluid was solved using the FDM to describe the plate bending and the BEM to describe the occurrence of fluid (FDM-BEM analysis).
The Gaussian uncertainty natural frequency is defined by its mean value and the specific range of its coefficient of variation, i.e., α(b) ∈ [0.00, 0.025]. The global response functions for first two natural frequencies were obtained in this case in the form of the following third-order polynomials: The results of the calculation with the expected values E(ω), coefficients of variation α(ω), skewness β(ω), and kurtosis κ(ω) of the first and second natural frequencies are presented in Figure 20.

Estimation of the Probabilistic Relative Entropy and the Safety Measure for Two Cantilever Plates Immersed in Fluid and Placed in One Plane along One Line
Probabilistic relative entropy H and the safety measure βH were estimated for the semi-analytical method (SAM) and stochastic perturbation technique (SPT), as well as the BEM-BEM, FEM-BEM, and FDM-BEM approaches and selected random parameters according to the relations (40) and (41), respectively.

Probabilistic Entropy for the Random Parameter Describing the Distance between Plates and the BEM-BEM Approach
The probabilistic relative entropy H estimated for the SAM and SPT solutions and the random distribution of the parameter a describing the distance between two plates for the first and second natural frequencies are presented in Figure 21.

Estimation of the Probabilistic Relative Entropy and the Safety Measure for Two Cantilever
Plates Immersed in Fluid and Placed in One Plane along One Line Probabilistic relative entropy H and the safety measure β H were estimated for the semi-analytical method (SAM) and stochastic perturbation technique (SPT), as well as the BEM-BEM, FEM-BEM, and FDM-BEM approaches and selected random parameters according to the relations (40) and (41), respectively.

Probabilistic Entropy for the Random Parameter Describing the Distance between Plates and the BEM-BEM Approach
The probabilistic relative entropy H estimated for the SAM and SPT solutions and the random distribution of the parameter a describing the distance between two plates for the first and second natural frequencies are presented in Figure 21.

Estimation of the Probabilistic Relative Entropy and the Safety Measure for Two Cantilever Plates Immersed in Fluid and Placed in One Plane along One Line
Probabilistic relative entropy H and the safety measure βH were estimated for the semi-analytical method (SAM) and stochastic perturbation technique (SPT), as well as the BEM-BEM, FEM-BEM, and FDM-BEM approaches and selected random parameters according to the relations (40) and (41), respectively.

Probabilistic Entropy for the Random Parameter Describing the Distance between Plates and the BEM-BEM Approach
The probabilistic relative entropy H estimated for the SAM and SPT solutions and the random distribution of the parameter a describing the distance between two plates for the first and second natural frequencies are presented in Figure 21.  The safety measure β H estimated for the random distribution of the same parameter a using the SAM and SPT solutions for the first and second natural frequencies are presented in Figure 22. The safety measure βH estimated for the random distribution of the same parameter a using the SAM and SPT solutions for the first and second natural frequencies are presented in Figure 22. It was observed that both methods returned practically the same results for the given variability range of the input uncertainty when studying the relative probabilistic entropy and its transformation to the reliability index. System safety exponentially decreases with an increase in the input CoV value; an analogous effect could be obtained while using the FORM index. A negligible effect of the given parameter uncertainty was observed in the limit values of both state functions for α(a) close to 0.25. The corresponding admissible value of the reliability index is close to 5, but equals about 40 in this example, which indicates a very small danger of resonance in this system.

Probabilistic Entropy for the Random Fluid Density and the FEM-BEM Approach
The probabilistic relative entropy H estimated for the SAM and SPT solutions with the random behavior of the fluid density for the first and second natural frequencies are presented in Figure 23. It was observed that both methods returned practically the same results for the given variability range of the input uncertainty when studying the relative probabilistic entropy and its transformation to the reliability index. System safety exponentially decreases with an increase in the input CoV value; an analogous effect could be obtained while using the FORM index. A negligible effect of the given parameter uncertainty was observed in the limit values of both state functions for α(a) close to 0.25. The corresponding admissible value of the reliability index is close to 5, but equals about 40 in this example, which indicates a very small danger of resonance in this system.

Probabilistic Entropy for the Random Fluid Density and the FEM-BEM Approach
The probabilistic relative entropy H estimated for the SAM and SPT solutions with the random behavior of the fluid density for the first and second natural frequencies are presented in Figure 23.
variability range of the input uncertainty when studying the relative probabilistic entropy and its transformation to the reliability index. System safety exponentially decreases with an increase in the input CoV value; an analogous effect could be obtained while using the FORM index. A negligible effect of the given parameter uncertainty was observed in the limit values of both state functions for α(a) close to 0.25. The corresponding admissible value of the reliability index is close to 5, but equals about 40 in this example, which indicates a very small danger of resonance in this system.

Probabilistic Entropy for the Random Fluid Density and the FEM-BEM Approach
The probabilistic relative entropy H estimated for the SAM and SPT solutions with the random behavior of the fluid density for the first and second natural frequencies are presented in Figure 23.  The safety measure β H evaluated for the random behavior of the fluid density by the SAM and SPT solutions for the first and second natural frequencies are presented in Figure 24. In this case study, the semi-analytical and perturbation-based methodologies returned almost the same values. It is clear from these results that the plate system safety is many times smaller when randomizing the fluid density. Even for very small values of this density coefficient of variation (close to 0, which is rather unusual in experimental fluid mechanics), the reliability index based on the Bhattacharyya relative entropy is below the well-known limits presented in the literature.

Probabilistic Entropy for the Random Fluid Density and the FDM-BEM Approach
The probabilistic relative entropy H estimated for the SAM and SPT solutions with the random behavior of the fluid density for the first and second natural frequencies are presented in Figure 25. In this case study, the semi-analytical and perturbation-based methodologies returned almost the same values. It is clear from these results that the plate system safety is many times smaller when randomizing the fluid density. Even for very small values of this density coefficient of variation (close to 0, which is rather unusual in experimental fluid mechanics), the reliability index based on the Bhattacharyya relative entropy is below the well-known limits presented in the literature.

Probabilistic Entropy for the Random Fluid Density and the FDM-BEM Approach
The probabilistic relative entropy H estimated for the SAM and SPT solutions with the random behavior of the fluid density for the first and second natural frequencies are presented in Figure 25.
low the well-known limits presented in the literature.

Probabilistic Entropy for the Random Fluid Density and the FDM-BEM Approach
The probabilistic relative entropy H estimated for the SAM and SPT solutions with the random behavior of the fluid density for the first and second natural frequencies are presented in Figure 25.  The safety measure β H estimated for the random behavior of the fluid density for the SAM and SPT solutions for the first and second natural frequencies are presented in Figure 26.

Concluding Remarks
The present paper presents the probabilistic approaches of Boundary Element, Finite Element, and Finite Difference Methods. Three probabilistic techniques were applied: the semi-analytical method (SAM), stochastic generalized perturbation technique (SPT), and Monte Carlo simulation (MCS).
The finite set of deterministic solutions necessary to proceed with probabilistic procedures was determined using the direct BEM approach with simplification of the formulation of the boundary and domain integral equations for isotropic and linear elastic thin (Kirchhoff-Love) plate theories. In addition, some cases were solved using the classic FEM and FDM approaches. The dynamic interaction between a plate and fluid was specified by the theory of the potential of the double layer for a fluid and the direct BEM approach.
From the viewpoint of the number of computational operations necessary to solve a single deterministic task, the optimal choice is the use of the FEM-BEM approach. In turn, the most regular arrangement of plate boundary elements and subsurfaces is only possible in terms of the BEM-BEM. In the other two cases (FEM-BEM and FDM-BEM), in order to maintain identical dimensions of the finite elements or the grid of points, the location of subsurfaces for liquids should always be adjusted so that the vertical component of the inertial force vector of the liquid coincides with the vertical component of the inertial force vector of the plate itself. Each of the approaches (BEM-BEM, FEM-BEM, and FDM-BEM) used requires the construction of a complete and generally asymmetric matrix H. This

Concluding Remarks
The present paper presents the probabilistic approaches of Boundary Element, Finite Element, and Finite Difference Methods. Three probabilistic techniques were applied: the semi-analytical method (SAM), stochastic generalized perturbation technique (SPT), and Monte Carlo simulation (MCS).
The finite set of deterministic solutions necessary to proceed with probabilistic procedures was determined using the direct BEM approach with simplification of the formulation of the boundary and domain integral equations for isotropic and linear elastic thin (Kirchhoff-Love) plate theories. In addition, some cases were solved using the classic FEM and FDM approaches. The dynamic interaction between a plate and fluid was specified by the theory of the potential of the double layer for a fluid and the direct BEM approach.
From the viewpoint of the number of computational operations necessary to solve a single deterministic task, the optimal choice is the use of the FEM-BEM approach. In turn, the most regular arrangement of plate boundary elements and subsurfaces is only possible in terms of the BEM-BEM. In the other two cases (FEM-BEM and FDM-BEM), in order to maintain identical dimensions of the finite elements or the grid of points, the location of subsurfaces for liquids should always be adjusted so that the vertical component of the inertial force vector of the liquid coincides with the vertical component of the inertial force vector of the plate itself. Each of the approaches (BEM-BEM, FEM-BEM, and FDM-BEM) used requires the construction of a complete and generally asymmetric matrix H. This matrix must then be inverted, which, with low accuracy, may cause numerical errors. Because of the large number of internal subsurfaces into which both surfaces of the plates are divided, it is time-consuming. From the viewpoint of the simplest formulation of the presented problem, the BEM-BEM approach seems to be the simplest and most useful. According to this approach, the division of the interior of the plates into subareas related to the liquid was the most regular and simple.
One can even be tempted to say that the joint implementation and application of SBEM, SFEM, and SFDM are not available in the literature in terms of any of the known stochastic methods; in this study, three approaches are used simultaneously. Third-degree polynomials were applied to fit the response signal curves for random parameters such as fluid density, plate placement, plate thickness, and collocation point location. It must be mentioned that the time consumption of both the semi-analytical and stochastic perturbation approaches are quasi-similar. The most time-consuming is the Monte Carlo simulation technique, which is attributable to the number of trials. The fluid density and boundary element collocation point placement (for BEM-BEM analysis) randomness seem to be negligible. From the viewpoint of the use of reliability measures expressed best by the coefficient β H , one can observe that their values drop sharply when the value of the argument, the uncertainty coefficient, increases.
The simplest of the random formulations-the SAM approach-allows the obtention of random moments in analytical forms, which is very convenient for further computational implementation in any academic and commercial BEM, FEM, or FDM programs.
The presented concept can also be successfully applied to the analysis of the dynamics of the control surfaces of ships or aircraft.